function Moment = elpm(Mean, Sigma, MAR, Order)
%ELPM Compute expected lower partial moments for a normal asset returns.
%	Given NUMSERIES asset returns with a vector of mean returns in a NUMSERIES
%	vector Mean, a vector of standard deviations of returns in a NUMSERIES
%	vector Sigma, a scalar minimum acceptable return MAR, and one or more
%	non-negative integer moment orders in a NUMORDERS vector Order, compute
%	expected lower partial moments relative to MAR for each asset in a
%	NUMORDERS x NUMSERIES matrix Moment.
%
%		elpm(Mean, Sigma);
%		elpm(Mean, Sigma, MAR)
%		elpm(Mean, Sigma, MAR, Order);
%		Moment = elpm(Mean, Sigma, MAR, Order);
%
% Inputs:
%	Mean - NUMSERIES vector with mean returns for a collection of NUMSERIES
%		assets.
%	Sigma - NUMSERIES vector with standard deviation of returns for a collection
%		of NUMSERIES assets.
%
% Optional Inputs:
%	MAR - Scalar minimum acceptable return (default MAR = 0). This is a cutoff
%		level of return such that all returns above MAR contribute nothing to
%		the lower partial moment.
%	Order - Either a scalar or an NUMORDERS-vector of non-negative integer
%		moment orders. If no order specified, default Order = 0, which is the
%		shortfall probability. This function will not work for negative or
%		non-integer orders.
%
% Outputs:
%	Moment - NUMORDERS x NUMSERIES matrix of expected lower partial moments with
%		NUMORDERS Orders and NUMSERIES series, i.e., each row contains
%		expected lower-partial moments for a given order.
%
% Notes:
%	To compute upper partial moments, just reverse the signs of both the input
%	Mean and MAR (do not reverse the signs of either Sigma or the output).
%
%	This function computes expected lower partial moments with the mean and
%	standard deviation of normally distributed asset returns. To compute sample
%	lower partial moments from asset returns, which has no distributional
%	assumptions, use lpm.
%
% See also lpm.

%	Copyright 1995-2006 The MathWorks, Inc.
%	$Revision: 1.1.6.2 $   $Date: 2006/06/16 20:09:40 $

% Step 1 - check arguments

if nargin < 2 || isempty(Mean) || isempty(Sigma)
	error('Finance:elpm:MissingInputArg', ...
		'Missing required input arguments Mean and Sigma.');
end

if isvector(Mean) && isa(Mean,'double')
	Mean = Mean(:);
	n = numel(Mean);
else
	error('Finance:elpm:InvalidInputArg', ...
 		'Invalid format for Mean. Must be a scalar or vector.');
end

if isvector(Sigma) && isa(Sigma,'double')
	Sigma = Sigma(:);
	ns = numel(Sigma);
	if ns ~= n
		error('Finance:elpm:InconsistentDims', ...
				'Extent of Mean (%d) and Sigma (%d) must be the same.',n,ns);
	end
else
	error('Finance:elpm:InvalidInputArg', ...
 		'Invalid format for Sigma. Must be a scalar or vector.');
end

if nargin < 3 || isempty(MAR)
	MAR = 0;
end

if ~isscalar(MAR) || ~isa(MAR, 'double')
	error('Finance:elpm:InvalidInputArg', ...
		'Invalid format for MAR. Must be a scalar.');
end

if nargin < 4 || isempty(Order)
	Order = 0;
else
	if isvector(Order) && isa(Order,'double')
		Order = Order(:);
	else
		error('Finance:elpm:InvalidInputArg', ...
			'Invalid format for Order. Must be a scalar or vector.');
	end
end

if any(floor(Order) ~= Order) || any(Order < 0)
	error('Finance:elpm:InvalidInputArg', ...
		'Invalid format for Order. Must be non-negative integers.');
end

p = numel(Order);

% Step 2 - compute expected lower partial moments

MaxOrder = max(Order);
Moment = zeros(p, n);

for i = 1:n
	AllMoments = elpmrange(Mean(i), Sigma(i), MAR, MaxOrder);
	for j = 1:p
		Moment(j,i) = AllMoments(Order(j)+1);
	end
end

% elpmrange - Recursion to compute all expected lower partial moments for all
%	orders 0, 1, ... , MaxOrder.

function Moments = elpmrange(Mean, Sigma, MAR, MaxOrder)

Moments = zeros(MaxOrder + 1,1);

if Sigma == 0
	if Mean <= MAR
		Moments = (MAR - Mean) .^ (0:MaxOrder);
	end
else
	Alpha = (Mean - MAR)/(Sigma*sqrt(2));

	Moments(1) = 0.5*erfc(Alpha);

	if (MaxOrder == 0)
		return
	end

	Moments(2) = Sigma*exp(-Alpha*Alpha)/sqrt(2*pi) - (Mean - MAR)*Moments(1);

	if (MaxOrder == 1)
		return
	end

	MM = Mean - MAR;
	SS = Sigma*Sigma;

	for i = 0:(MaxOrder - 2);
		Moments(i + 3) = (i + 1)*SS*Moments(i + 1) - MM*Moments(i + 2);
	end
end
